################################################################################
################################ PLOT ##########################################
################################################################################
##########Load results for a simulation with 200 observations per group.########
k=4
numobs=200

fname <- paste("C:\\DASEV\\S1\\simresult_N_",numobs,"_",k,".Rdata",sep="")
fname1 <- paste("C:\\DASEV\\S1\\otherresult_N_",numobs,"_",k,".Rdata",sep="")
fnamep <- paste("C:\\DASEV\\S1\\paradata_N_",numobs,"_",k,".Rdata",sep="")
fname.data <- paste("C:\\DASEV\\S1\\simdata_N_",numobs,"_",k,".Rdata",sep="")

load(file=fname.data) # data2
datatmp <- data2
count <- rowSums(datatmp==0)
rawpall  <- count/(2*numobs)

load(file=fname) #sim.result
est <- sim.result$estimates
p <-  sim.result$pvalue_both
p.m <-  sim.result$pvalue_mean
p.p <-  sim.result$pvalue_zero
post.mucontrol.f <- est[,1]
post.mucase.f <- est[,1]+est[,2]
post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
post.var <- est[,5]^2

load(file=fname1) #other.result
t.p <- other.result[,1]
t.var <- other.result[,2]
t.mucontrol <- other.result[,3]
t.mucase <- other.result[,4]
t.pcontrol <- other.result[,5]
t.pcase <- other.result[,6]
t.p.m <- other.result[,7]
t.p.p <- other.result[,13]

load(file=fnamep) #parameters
sim.lod <- paradata[,1]
sim.mucontrol <- paradata[,2]
sim.pcontrol <- paradata[,3]
sim.sd <- paradata[,4]
sim.mucase <- paradata[,5]
sim.pcase <- paradata[,6]
sim.geneid <- paradata[,7]
sim.geneid1 <- paradata[,8]
numint <- paradata[,9]

sim.var <- sim.sd^2
q<-p.adjust(p,method='BH')
qmean<-p.adjust(p.m,method='BH')
t.q<-p.adjust(t.p,method='BH')
t.qmean<-p.adjust(t.p.m,method="BH")

################################Figure 1#########################################
par(mfrow=c(1,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))
rank <- c(1:length(sim.geneid))

sim.sort.result1 <- cbind(sim.geneid,qmean,post.mucase.f,post.mucontrol.f,post.pcase.f,post.pcontrol.f,post.var,numint,p.m,sim.pcontrol)
sim.sort.result1 <-sim.sort.result1[order(p.m),]
sim.sort.result1 <- cbind(sim.sort.result1,rank)

t.sort.result1 <- cbind(sim.geneid,t.qmean,t.mucase,t.mucontrol,t.pcase,t.pcontrol,t.var,numint,t.p.m,sim.pcontrol)
t.sort.result1 <-t.sort.result1[order(t.p.m),]
t.sort.result1 <- cbind(t.sort.result1,rank)

TOPN = 150
col1=rgb(0.1,0.3,0.5,1)
col3=rgb(0.8,0.3,0,1)

resultsub <- t.sort.result1[1:TOPN,]
sim.logfc <- resultsub[,3]-resultsub[,4]
shapint<-ifelse(resultsub[,1]==1,1,0)
realNpoints <- shapint == 0
realPpoints <- shapint == 1
zeroint <- resultsub[,10]>0.9
plot(1, type="n", xlab="Log Fold Change", ylab="Variance", xlim=c(-4,4), ylim=c(0, max(resultsub[,7]+1)),main="a",cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(sim.logfc[realPpoints], resultsub[,7][realPpoints],col=col1,pch=2,cex =1)
points(sim.logfc[realNpoints], resultsub[,7][realNpoints],col=col3,pch=1)
points(sim.logfc[realPpoints&zeroint], resultsub[,7][realPpoints&zeroint],col=col1,pch=17,cex =1)
points(sim.logfc[realNpoints&zeroint], resultsub[,7][realNpoints&zeroint],col=col3,pch=16)
legend("topright", legend=c("FP PMV <= 90%","TP PMV <= 90%","FP PMV > 90%","TP PMV > 90%"),border="white", col=c(c(col3,col1),c(col3,col1)),pch=c(c(1,2),c(16,17)),cex=1.2)

resultsub <- sim.sort.result1[1:TOPN,]
sim.logfc <- resultsub[,3]-resultsub[,4]
shapint<-ifelse(resultsub[,1]==1,1,0)
realNpoints <- shapint == 0
realPpoints <- shapint == 1
zeroint <- resultsub[,10]>0.9
plot(1, type="n", xlab="Log Fold Change", ylab="Variance", xlim=c(-4,4), ylim=c(0, max(resultsub[,7]+1)),main="b",cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(sim.logfc[realPpoints], resultsub[,7][realPpoints],col=col1,pch=2,cex =1)
points(sim.logfc[realNpoints], resultsub[,7][realNpoints],col=col3,pch=1)
points(sim.logfc[realPpoints&zeroint], resultsub[,7][realPpoints&zeroint],col=col1,pch=17,cex =1)
points(sim.logfc[realNpoints&zeroint], resultsub[,7][realNpoints&zeroint],col=col3,pch=16)

################################Figure 2########################################
par(mfrow=c(2,2))
par(oma=c(2,2,0,0))
par(mar=c(3,3,2,1))

pind <- rawpall>0.9
pind2 <- abs(post.var-sim.var)<0.5
pind3 <- abs(t.var-sim.var)<0.6
pind4 <- abs(post.var-sim.var)<0.2
pind5 <- abs(t.var-sim.var)<0.2
cexsize=0.6
col1=rgb(0.1,0.3,0.5,0.8)
col2=rgb(0.8,0.3,0,0.1)
col3=rgb(0.8,0.3,0,0.6)
col4=rgb(0.8,0.3,0,0.3)
col5=rgb(0.5,0.3,0,0.1)
par(mfrow=c(2,2))
par(oma=c(2,2,0,0))
par(mar=c(3,3,2,1))
plot(t.var[pind==FALSE],sim.var[pind==FALSE],xlim= c(0,22),xlab="", ylab="",main = "a",cex.axis=1.4,cex.main=2,pch=16,col=col1,cex =cexsize)
points(t.var[pind==T],sim.var[pind==T],pch=16,col=col2,cex =cexsize)
points(t.var[pind==T&pind3==FALSE],sim.var[pind==T&pind3==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
plot(post.var[pind==FALSE],sim.var[pind==FALSE],xlim= c(0,22),xlab="", ylab="",main = "b",cex.axis=1.4,cex.lab=2,cex.main=2,pch=16,col=col1,cex =cexsize)
points(post.var[pind==T],sim.var[pind==T],pch=16,col=col5,cex =cexsize)
points(post.var[pind==T&pind2==FALSE],sim.var[pind==T&pind2==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
plot(t.var[pind==FALSE],sim.var[pind==FALSE],xlim=c(0,2), ylim=c(0,2),xlab="", ylab="",main = "c",cex.axis=1.4,cex.lab=2,cex.main=2,pch=16,col=col1,cex =cexsize)
points(t.var[pind==T],sim.var[pind==T],pch=16,col=col4,cex =cexsize)
points(t.var[pind==T&pind5==FALSE],sim.var[pind==T&pind5==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
plot(post.var[pind==FALSE],sim.var[pind==FALSE],xlim=c(0,2), ylim=c(0,2),xlab="", ylab="",main = "d",cex.axis=1.4,cex.main=2,pch=16,col=col1,cex =cexsize)
points(post.var[pind==T],sim.var[pind==T],pch=16,col=col4,cex =cexsize)
points(post.var[pind==T&pind4==FALSE],sim.var[pind==T&pind4==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
legend("bottomright", legend=c(paste("PMV > 90%"),paste("PMV <= 90%")),border="white", col=c(col3,col1),pch=c(16,16),cex=1.3)

mtext("True Variance", side=2, line=0, cex=2,outer=TRUE)
mtext("Estimated Variance", side=1, line=0, cex=2,outer=TRUE)

################################Figure 3########################################
par(mfrow=c(2,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))

cexsize=0.6
col1=rgb(0.1,0.3,0.5,0.8)
col2=rgb(0.8,0.3,0,0.1)
col3=rgb(0.8,0.3,0,0.6)
col4=rgb(0.8,0.3,0,0.3)
col5=rgb(0.5,0.3,0,0.1)
par(mfrow=c(2,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))
pind <- rawpall>0.9
pind2<- abs(t.mucontrol-sim.mucontrol)>0.5
pind3<- abs(post.mucontrol.f-sim.mucontrol)>0.5
pind4<- abs(t.pcontrol-sim.pcontrol)>0.1
pind5<- abs(post.pcontrol.f-sim.pcontrol)>0.1
plot(1, type="n", xlab="Estimated Mean", ylab="True Mean",main = "a",xlim=c(0,10), ylim=c(0,10),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(t.mucontrol[pind==FALSE], sim.mucontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(t.mucontrol[pind], sim.mucontrol[pind],col=col2,pch=16,cex=cexsize)
points(t.mucontrol[pind&pind2], sim.mucontrol[pind&pind2],col=col3,pch=16,cex=cexsize)
plot(1, type="n", xlab="Estimated Mean", ylab="True Mean",main = "b",xlim=c(0,10), ylim=c(0,10),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(post.mucontrol.f[pind==FALSE], sim.mucontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(post.mucontrol.f[pind], sim.mucontrol[pind],col=col2,pch=16,cex=cexsize)
points(post.mucontrol.f[pind&pind3], sim.mucontrol[pind&pind3],col=col3,pch=16,cex=cexsize)
#legend("bottomright", legend=c(paste("PMV > 90%"),paste("PMV <= 90%")),border="white", col=c(col3,col1),pch=c(16,16),cex=1.3)

plot(1, type="n", xlab="Estimated BPMV Proprotion", ylab="True BPMV Proportion",main = "c",xlim=c(0,1), ylim=c(0,1),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(t.pcontrol[pind==FALSE], sim.pcontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(t.pcontrol[pind], sim.pcontrol[pind],col=col2,pch=16,cex=cexsize)
points(t.pcontrol[pind&pind4], sim.pcontrol[pind&pind4],col=col3,pch=16,cex=cexsize)
plot(1, type="n", xlab="Estimated BPMV Proprotion", ylab="True BPMV Proportion",main = "d",xlim=c(0,1), ylim=c(0,1),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(post.pcontrol.f[pind==FALSE], sim.pcontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(post.pcontrol.f[pind], sim.pcontrol[pind],col=col2,pch=16,cex=cexsize)
points(post.pcontrol.f[pind&pind5], sim.pcontrol[pind&pind5],col=col3,pch=16,cex=cexsize)
legend("bottomright", legend=c(paste("PMV > 90%"),paste("PMV <= 90%")),border="white", col=c(col3,col1),pch=c(16,16),cex=1.3)

################################Figure 4########################################
#Average Reported FDR vs True FDR over 100 simulation
par(mfrow=c(2,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))

numobs=100
DASEV.TRP <-c()
DASEV.TRP.M <-c()
LIM.TRP <-c()
LIM.TRP.M <-c()
DASEV.FDR <-c()
DASEV.FDR.M <-c()
LIM.FDR <-c()
LIM.FDR.M <-c()
a1_all <- c()
b1_all <- c()
c1_all <- c()
d1_all <- c()
a2_all <- c()
b2_all <- c()
c2_all <- c()
d2_all <- c()
a3_all <- c()
b3_all <- c()
c3_all <- c()
d3_all <- c()
for (k in 1:100){
  
  fname <- paste("C:\\DASEV\\S1\\simresult_N_",numobs,"_",k,".Rdata",sep="")
  fname1 <- paste("C:\\DASEV\\S1\\otherresult_N_",numobs,"_",k,".Rdata",sep="")
  fnamep <- paste("C:\\DASEV\\S1\\paradata_N_",numobs,"_",k,".Rdata",sep="")
  fname.data <- paste("C:\\DASEV\\S1\\simdata_N_",numobs,"_",k,".Rdata",sep="")
  
  load(file=fname.data) # data2
  datatmp <- data2
  count <- rowSums(datatmp==0)
  rawpall  <- count/(2*numobs)
  
  load(file=fname) #sim.result
  est <- sim.result$estimates
  p <-  sim.result$pvalue_both
  p.m <-  sim.result$pvalue_mean
  p.p <-  sim.result$pvalue_zero
  post.mucontrol.f <- est[,1]
  post.mucase.f <- est[,1]+est[,2]
  post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
  post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
  post.var <- est[,5]^2
  
  load(file=fname1) #other.result
  t.p <- other.result[,1]
  t.var <- other.result[,2]
  t.mucontrol <- other.result[,3]
  t.mucase <- other.result[,4]
  t.pcontrol <- other.result[,5]
  t.pcase <- other.result[,6]
  t.p.m <- other.result[,7]
  t.p.p <- other.result[,13]
  
  load(file=fnamep) #parameters
  sim.lod <- paradata[,1]
  sim.mucontrol <- paradata[,2]
  sim.pcontrol <- paradata[,3]
  sim.sd <- paradata[,4]
  sim.mucase <- paradata[,5]
  sim.pcase <- paradata[,6]
  sim.geneid <- paradata[,7]
  sim.geneid1 <- paradata[,8]
  numint <- paradata[,9]
  
  sim.var <- sim.sd^2
  q<-p.adjust(p,method='BH')
  qmean<-p.adjust(p.m,method='BH')
  t.q<-p.adjust(t.p,method='BH')
  t.qmean<-p.adjust(t.p.m,method="BH")
  
  rank <- c(1:length(sim.geneid))
  
  a1_ <- sum(qmean<=0.01 & sim.geneid==1)
  b1_ <- sum(qmean<=0.01 & sim.geneid==0)
  c1_ <- sum(t.qmean<=0.01 & sim.geneid==1)
  d1_ <- sum(t.qmean<=0.01 & sim.geneid==0)
  
  a2_ <- sum(qmean<=0.05 & sim.geneid==1)
  b2_ <- sum(qmean<=0.05 & sim.geneid==0)
  c2_ <- sum(t.qmean<=0.05 & sim.geneid==1)
  d2_ <- sum(t.qmean<=0.05 & sim.geneid==0)
  
  a3_ <- sum(qmean<=0.1 & sim.geneid==1)
  b3_ <- sum(qmean<=0.1 & sim.geneid==0)
  c3_ <- sum(t.qmean<=0.1 & sim.geneid==1)
  d3_ <- sum(t.qmean<=0.1 & sim.geneid==0)
  
  a1_all <- c(a1_all,a1_)
  b1_all <- c(b1_all,b1_)
  c1_all <- c(c1_all,c1_)
  d1_all <- c(d1_all,d1_)
  a2_all <- c(a2_all,a2_)
  b2_all <- c(b2_all,b2_)
  c2_all <- c(c2_all,c2_)
  d2_all <- c(d2_all,d2_)
  a3_all <- c(a3_all,a3_)
  b3_all <- c(b3_all,b3_)
  c3_all <- c(c3_all,c3_)
  d3_all <- c(d3_all,d3_)
  
  sim.sort.result <- cbind(sim.geneid,q,post.mucase.f,numint,p)
  sim.sort.result <-sim.sort.result[order(p),]
  sim.sort.result <- cbind(sim.sort.result,rank)
  sim.sort.result1 <- cbind(sim.geneid,qmean,post.mucase.f,numint,p.m)
  sim.sort.result1 <-sim.sort.result1[order(p.m),]
  sim.sort.result1 <- cbind(sim.sort.result1,rank)
  t.sort.result <- cbind(sim.geneid,t.q,t.mucase,numint,t.p)
  t.sort.result <-t.sort.result[order(t.p),]
  t.sort.result <- cbind(t.sort.result,rank)
  t.sort.result1 <- cbind(sim.geneid,t.qmean,t.mucase,numint,t.p.m)
  t.sort.result1 <-t.sort.result1[order(t.p.m),]
  t.sort.result1 <- cbind(t.sort.result1,rank)
  
  #tOP RANKED FEATURES AND TRUE POSITIVE RATE
  TOPn <- seq(5,2000,5)
  TPR <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR <- c(TPR,tpr)
    TPR[is.na(TPR)] <- 0
  }
  DASEV.TRP <-cbind(DASEV.TRP,TPR)
  
  TPR.M <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR.M <- c(TPR.M,tpr)
    TPR.M[is.na(TPR.M)] <- 0
  }
  DASEV.TRP.M <-cbind(DASEV.TRP.M,TPR.M)
  
  T.TPR <- c()
  for (i in TOPn){
    resultsub <- t.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR <- c(T.TPR,tpr)
    T.TPR[is.na(T.TPR)] <- 0
  }
  LIM.TRP <- cbind(LIM.TRP,T.TPR)
  
  T.TPR.M <- c()
  for (i in TOPn){
    resultsub <- t.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR.M <- c(T.TPR.M,tpr)
    T.TPR.M[is.na(T.TPR.M)] <- 0
  }
  LIM.TRP.M <- cbind(LIM.TRP.M,T.TPR.M)
  
  cutoffs <- seq(0,1,0.005)
  
  FDR <- c()
  FDR.M <- c()
  T.FDR <- c()
  T.FDR.M <- c()
  for (j in cutoffs){
    DE=sum(q<j & sim.geneid==1)
    TOT=sum(q<j)
    if (TOT==0) fdr=0 else fdr = 1- DE/TOT
    FDR <- c(FDR,fdr)
    
    DE.M=sum(qmean<j & sim.geneid==1)
    TOT.M=sum(qmean<j)
    if (TOT.M==0) fdr=0 else fdr = 1- DE.M/TOT.M
    FDR.M <- c(FDR.M,fdr)
    
    T.DE=sum(t.q<j & sim.geneid==1)
    T.TOT=sum(t.q<j)
    if (T.TOT==0) fdr=0 else fdr = 1- T.DE/T.TOT
    T.FDR <- c(T.FDR,fdr)
    
    T.DE.M=sum(t.qmean<j & sim.geneid==1)
    T.TOT.M=sum(t.qmean<j)
    if (T.TOT.M==0) fdr=0 else fdr = 1- T.DE.M/T.TOT.M
    T.FDR.M <- c(T.FDR.M,fdr)
  }
  DASEV.FDR <-cbind(DASEV.FDR,FDR)
  DASEV.FDR.M <-cbind(DASEV.FDR.M,FDR.M)
  LIM.FDR <-cbind(LIM.FDR,T.FDR)
  LIM.FDR.M <-cbind(LIM.FDR.M,T.FDR.M)
}

mean.DASEV.TRP_100 <-rowMeans(DASEV.TRP)
mean.DASEV.TRP.M_100  <-rowMeans(DASEV.TRP.M)
mean.LIM.TRP_100  <-rowMeans(LIM.TRP)
mean.LIM.TRP.M_100  <-rowMeans(LIM.TRP.M)
mean.DASEV.FDR_100  <-rowMeans(DASEV.FDR)
mean.DASEV.FDR.M_100  <-rowMeans(DASEV.FDR.M)
mean.LIM.FDR_100  <-rowMeans(LIM.FDR)
mean.LIM.FDR.M_100  <-rowMeans(LIM.FDR.M)
a1<- round(mean(a1_all))
b1<- round(mean(b1_all))
c1<- round(mean(c1_all))
d1<- round(mean(d1_all))
a2<- round(mean(a2_all))
b2<- round(mean(b2_all))
c2<- round(mean(c2_all))
d2<- round(mean(d2_all))
a3<- round(mean(a3_all))
b3<- round(mean(b3_all))
c3<- round(mean(c3_all))
d3<- round(mean(d3_all))

v1_100  <- c(a1,a2,a3)
v2_100  <- c(-b1,-b2,-b3)
v3_100  <- c(c1,c2,c3)
v4_100  <- c(-d1,-d2,-d3)
#  cutoffs2<- c(0.01,0.05,0.1)
tfdrd1 <- b1_all/(b1_all+a1_all)
tfdrd2 <- b2_all/(b2_all+a2_all)
tfdrd3 <- b3_all/(b3_all+a3_all)
tfdrt1 <- d1_all/(d1_all+c1_all)
tfdrt2 <- d2_all/(d2_all+c2_all)
tfdrt3 <- d3_all/(d3_all+c3_all)

tprrate_100_2  <- c(mean(tfdrd1[!is.na(tfdrd1)]),mean(tfdrt1[!is.na(tfdrt1)]),
                    mean(tfdrd2[!is.na(tfdrd2)]),mean(tfdrt2[!is.na(tfdrt2)]),
                    mean(tfdrd3[!is.na(tfdrd3)]),mean(tfdrt3[!is.na(tfdrd3)]))
rate_100  <- paste(round(100*tprrate_100_2 , 2), "%", sep="")

numobs=200
DASEV.TRP <-c()
DASEV.TRP.M <-c()
LIM.TRP <-c()
LIM.TRP.M <-c()
DASEV.FDR <-c()
DASEV.FDR.M <-c()
LIM.FDR <-c()
LIM.FDR.M <-c()
a1_all <- c()
b1_all <- c()
c1_all <- c()
d1_all <- c()
a2_all <- c()
b2_all <- c()
c2_all <- c()
d2_all <- c()
a3_all <- c()
b3_all <- c()
c3_all <- c()
d3_all <- c()
for (k in 1:100){
  #School PC
  
  fname <- paste("C:\\DASEV\\S1\\simresult_N_",numobs,"_",k,".Rdata",sep="")
  fname1 <- paste("C:\\DASEV\\S1\\otherresult_N_",numobs,"_",k,".Rdata",sep="")
  fnamep <- paste("C:\\DASEV\\S1\\paradata_N_",numobs,"_",k,".Rdata",sep="")
  fname.data <- paste("C:\\DASEV\\S1\\simdata_N_",numobs,"_",k,".Rdata",sep="")
  
  load(file=fname.data) # data2
  datatmp <- data2
  count <- rowSums(datatmp==0)
  rawpall  <- count/(2*numobs)
  
  load(file=fname) #sim.result
  est <- sim.result$estimates
  p <-  sim.result$pvalue_both
  p.m <-  sim.result$pvalue_mean
  p.p <-  sim.result$pvalue_zero
  post.mucontrol.f <- est[,1]
  post.mucase.f <- est[,1]+est[,2]
  post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
  post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
  post.var <- est[,5]^2
  
  load(file=fname1) #other.result
  t.p <- other.result[,1]
  t.var <- other.result[,2]
  t.mucontrol <- other.result[,3]
  t.mucase <- other.result[,4]
  t.pcontrol <- other.result[,5]
  t.pcase <- other.result[,6]
  t.p.m <- other.result[,7]
  t.p.p <- other.result[,13]
  
  load(file=fnamep) #parameters
  sim.lod <- paradata[,1]
  sim.mucontrol <- paradata[,2]
  sim.pcontrol <- paradata[,3]
  sim.sd <- paradata[,4]
  sim.mucase <- paradata[,5]
  sim.pcase <- paradata[,6]
  sim.geneid <- paradata[,7]
  sim.geneid1 <- paradata[,8]
  numint <- paradata[,9]
  
  sim.var <- sim.sd^2
  q<-p.adjust(p,method='BH')
  qmean<-p.adjust(p.m,method='BH')
  t.q<-p.adjust(t.p,method='BH')
  t.qmean<-p.adjust(t.p.m,method="BH")
  
  rank <- c(1:length(sim.geneid))
  
  a1_ <- sum(qmean<=0.01 & sim.geneid==1)
  b1_ <- sum(qmean<=0.01 & sim.geneid==0)
  c1_ <- sum(t.qmean<=0.01 & sim.geneid==1)
  d1_ <- sum(t.qmean<=0.01 & sim.geneid==0)
  
  a2_ <- sum(qmean<=0.05 & sim.geneid==1)
  b2_ <- sum(qmean<=0.05 & sim.geneid==0)
  c2_ <- sum(t.qmean<=0.05 & sim.geneid==1)
  d2_ <- sum(t.qmean<=0.05 & sim.geneid==0)
  
  a3_ <- sum(qmean<=0.1 & sim.geneid==1)
  b3_ <- sum(qmean<=0.1 & sim.geneid==0)
  c3_ <- sum(t.qmean<=0.1 & sim.geneid==1)
  d3_ <- sum(t.qmean<=0.1 & sim.geneid==0)
  
  a1_all <- c(a1_all,a1_)
  b1_all <- c(b1_all,b1_)
  c1_all <- c(c1_all,c1_)
  d1_all <- c(d1_all,d1_)
  a2_all <- c(a2_all,a2_)
  b2_all <- c(b2_all,b2_)
  c2_all <- c(c2_all,c2_)
  d2_all <- c(d2_all,d2_)
  a3_all <- c(a3_all,a3_)
  b3_all <- c(b3_all,b3_)
  c3_all <- c(c3_all,c3_)
  d3_all <- c(d3_all,d3_)
  
  
  sim.sort.result <- cbind(sim.geneid,q,post.mucase.f,numint,p)
  sim.sort.result <-sim.sort.result[order(p),]
  sim.sort.result <- cbind(sim.sort.result,rank)
  sim.sort.result1 <- cbind(sim.geneid,qmean,post.mucase.f,numint,p.m)
  sim.sort.result1 <-sim.sort.result1[order(p.m),]
  sim.sort.result1 <- cbind(sim.sort.result1,rank)
  t.sort.result <- cbind(sim.geneid,t.q,t.mucase,numint,t.p)
  t.sort.result <-t.sort.result[order(t.p),]
  t.sort.result <- cbind(t.sort.result,rank)
  t.sort.result1 <- cbind(sim.geneid,t.qmean,t.mucase,numint,t.p.m)
  t.sort.result1 <-t.sort.result1[order(t.p.m),]
  t.sort.result1 <- cbind(t.sort.result1,rank)
  
  #tOP RANKED FEATURES AND TRUE POSITIVE RATE
  TOPn <- seq(5,2000,5)
  TPR <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR <- c(TPR,tpr)
    TPR[is.na(TPR)] <- 0
  }
  DASEV.TRP <-cbind(DASEV.TRP,TPR)
  
  TPR.M <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR.M <- c(TPR.M,tpr)
    TPR.M[is.na(TPR.M)] <- 0
  }
  DASEV.TRP.M <-cbind(DASEV.TRP.M,TPR.M)
  
  T.TPR <- c()
  for (i in TOPn){
    resultsub <- t.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR <- c(T.TPR,tpr)
    T.TPR[is.na(T.TPR)] <- 0
  }
  LIM.TRP <- cbind(LIM.TRP,T.TPR)
  
  T.TPR.M <- c()
  for (i in TOPn){
    resultsub <- t.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR.M <- c(T.TPR.M,tpr)
    T.TPR.M[is.na(T.TPR.M)] <- 0
  }
  LIM.TRP.M <- cbind(LIM.TRP.M,T.TPR.M)
  
  
  cutoffs <- seq(0,1,0.005)
  
  FDR <- c()
  FDR.M <- c()
  T.FDR <- c()
  T.FDR.M <- c()
  for (j in cutoffs){
    DE=sum(q<j & sim.geneid==1)
    TOT=sum(q<j)
    if (TOT==0) fdr=0 else fdr = 1- DE/TOT
    FDR <- c(FDR,fdr)
    
    DE.M=sum(qmean<j & sim.geneid==1)
    TOT.M=sum(qmean<j)
    if (TOT.M==0) fdr=0 else fdr = 1- DE.M/TOT.M
    FDR.M <- c(FDR.M,fdr)
    
    T.DE=sum(t.q<j & sim.geneid==1)
    T.TOT=sum(t.q<j)
    if (T.TOT==0) fdr=0 else fdr = 1- T.DE/T.TOT
    T.FDR <- c(T.FDR,fdr)
    
    T.DE.M=sum(t.qmean<j & sim.geneid==1)
    T.TOT.M=sum(t.qmean<j)
    if (T.TOT.M==0) fdr=0 else fdr = 1- T.DE.M/T.TOT.M
    T.FDR.M <- c(T.FDR.M,fdr)
  }
  DASEV.FDR <-cbind(DASEV.FDR,FDR)
  DASEV.FDR.M <-cbind(DASEV.FDR.M,FDR.M)
  LIM.FDR <-cbind(LIM.FDR,T.FDR)
  LIM.FDR.M <-cbind(LIM.FDR.M,T.FDR.M)
}

mean.DASEV.TRP <-rowMeans(DASEV.TRP)
mean.DASEV.TRP.M <-rowMeans(DASEV.TRP.M)
mean.LIM.TRP <-rowMeans(LIM.TRP)
mean.LIM.TRP.M <-rowMeans(LIM.TRP.M)
mean.DASEV.FDR <-rowMeans(DASEV.FDR)
mean.DASEV.FDR.M <-rowMeans(DASEV.FDR.M)
mean.LIM.FDR <-rowMeans(LIM.FDR)
mean.LIM.FDR.M <-rowMeans(LIM.FDR.M)
a1<- round(mean(a1_all))
b1<- round(mean(b1_all))
c1<- round(mean(c1_all))
d1<- round(mean(d1_all))
a2<- round(mean(a2_all))
b2<- round(mean(b2_all))
c2<- round(mean(c2_all))
d2<- round(mean(d2_all))
a3<- round(mean(a3_all))
b3<- round(mean(b3_all))
c3<- round(mean(c3_all))
d3<- round(mean(d3_all))

plot(TOPn,mean.DASEV.TRP.M_100 ,type="l",col="darkorange2", xlim=c(0,300),ylim=c(0.5,1), xlab="Top Ranked Features", ylab="True Positive Rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(TOPn,mean.LIM.TRP.M_100 ,type="l",col="royalblue4", lwd=2)
legend("bottomleft", legend=c("DASEV","TLK"),border="white", col=c("darkorange2","royalblue4"),lty=c(1,1),cex=1.2)

plot(TOPn,mean.DASEV.TRP.M,type="l",col="darkorange2", xlim=c(0,300),ylim=c(0.5,1), xlab="Top Ranked Features", ylab="True Positive Rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(TOPn,mean.LIM.TRP.M,type="l",col="royalblue4", lwd=2)
legend("bottomleft", legend=c("DASEV","TLK"),border="white", col=c("darkorange2","royalblue4"),lty=c(1,1),cex=1.2)

barplot(matrix(c(v1_100 ,v3_100 ),nrow=2, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4'),ylim=c(-150,200),xlab="Reported FDR", ylab = "Counts",main = "c",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(v2_100 ,v4_100 ),nrow=2, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2'), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,4.5,5.5,7.5,8.5), y=c(90,110,120,140,150,170), labels=rate_100 ,cex=1.2)
legend("bottomleft", legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2"),cex=1.2,pch=c(15,15,15,15))

v1 <- c(a1,a2,a3)
v2 <- c(-b1,-b2,-b3)
v3 <- c(c1,c2,c3)
v4 <- c(-d1,-d2,-d3)
cutoffs2<- c(0.01,0.05,0.1)
tfdrd1 <- b1_all/(b1_all+a1_all)
tfdrd2 <- b2_all/(b2_all+a2_all)
tfdrd3 <- b3_all/(b3_all+a3_all)
tfdrt1 <- d1_all/(d1_all+c1_all)
tfdrt2 <- d2_all/(d2_all+c2_all)
tfdrt3 <- d3_all/(d3_all+c3_all)

tprrate_2  <- c(mean(tfdrd1[!is.na(tfdrd1)]),mean(tfdrt1[!is.na(tfdrt1)]),
                mean(tfdrd2[!is.na(tfdrd2)]),mean(tfdrt2[!is.na(tfdrt2)]),
                mean(tfdrd3[!is.na(tfdrd3)]),mean(tfdrt3[!is.na(tfdrd3)]))
rate <- paste(round(100*tprrate_2, 2), "%", sep="")
barplot(matrix(c(v1,v3),nrow=2, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4'),ylim=c(-200,350),xlab="Reported FDR", ylab = "Counts",main = "d",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(v2,v4),nrow=2, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2'), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,4.5,5.5,7.5,8.5), y=(c(150,180,190,220,230,260)+20), labels=rate,cex=1.2)
legend("bottomleft", legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2"),cex=1.2,pch=c(15,15,15,15))


#########Load results for a simulation with 100 observations per group.#########
k=2
numobs=100

fname <- paste("C:\\DASEV\\S1\\simresult_N_",numobs,"_",k,".Rdata",sep="")
fname1 <- paste("C:\\DASEV\\S1\\otherresult_N_",numobs,"_",k,".Rdata",sep="")
fnamep <- paste("C:\\DASEV\\S1\\paradata_N_",numobs,"_",k,".Rdata",sep="")
fname.data <- paste("C:\\DASEV\\S1\\simdata_N_",numobs,"_",k,".Rdata",sep="")

load(file=fname.data) # data2
datatmp <- data2
count <- rowSums(datatmp==0)
rawpall  <- count/(2*numobs)

load(file=fname) #sim.result
est <- sim.result$estimates
p <-  sim.result$pvalue_both
p.m <-  sim.result$pvalue_mean
p.p <-  sim.result$pvalue_zero
post.mucontrol.f <- est[,1]
post.mucase.f <- est[,1]+est[,2]
post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
post.var <- est[,5]^2

load(file=fname1) #other.result
t.p <- other.result[,1]
t.var <- other.result[,2]
t.mucontrol <- other.result[,3]
t.mucase <- other.result[,4]
t.pcontrol <- other.result[,5]
t.pcase <- other.result[,6]
t.p.m <- other.result[,7]
t.p.p <- other.result[,13]

load(file=fnamep) #parameters
sim.lod <- paradata[,1]
sim.mucontrol <- paradata[,2]
sim.pcontrol <- paradata[,3]
sim.sd <- paradata[,4]
sim.mucase <- paradata[,5]
sim.pcase <- paradata[,6]
sim.geneid <- paradata[,7]
sim.geneid1 <- paradata[,8]
numint <- paradata[,9]

sim.var <- sim.sd^2
q<-p.adjust(p,method='BH')
qmean<-p.adjust(p.m,method='BH')
t.q<-p.adjust(t.p,method='BH')
t.qmean<-p.adjust(t.p.m,method="BH")

################################Figure S2#########################################
par(mfrow=c(1,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))
rank <- c(1:length(sim.geneid))

sim.sort.result1 <- cbind(sim.geneid,qmean,post.mucase.f,post.mucontrol.f,post.pcase.f,post.pcontrol.f,post.var,numint,p.m,sim.pcontrol)
sim.sort.result1 <-sim.sort.result1[order(p.m),]
sim.sort.result1 <- cbind(sim.sort.result1,rank)

t.sort.result1 <- cbind(sim.geneid,t.qmean,t.mucase,t.mucontrol,t.pcase,t.pcontrol,t.var,numint,t.p.m,sim.pcontrol)
t.sort.result1 <-t.sort.result1[order(t.p.m),]
t.sort.result1 <- cbind(t.sort.result1,rank)

TOPN = 150
col1=rgb(0.1,0.3,0.5,1)
col3=rgb(0.8,0.3,0,1)

resultsub <- t.sort.result1[1:TOPN,]
sim.logfc <- resultsub[,3]-resultsub[,4]
shapint<-ifelse(resultsub[,1]==1,1,0)
realNpoints <- shapint == 0
realPpoints <- shapint == 1
zeroint <- resultsub[,10]>0.9
plot(1, type="n", xlab="Log Fold Change", ylab="Variance", xlim=c(-4,4), ylim=c(0, 4),main="a",cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(sim.logfc[realPpoints], resultsub[,7][realPpoints],col=col1,pch=2,cex =1)
points(sim.logfc[realNpoints], resultsub[,7][realNpoints],col=col3,pch=1)
points(sim.logfc[realPpoints&zeroint], resultsub[,7][realPpoints&zeroint],col=col1,pch=17,cex =1)
points(sim.logfc[realNpoints&zeroint], resultsub[,7][realNpoints&zeroint],col=col3,pch=16)
#legend("topright", legend=c(paste("FP =",sum(realNpoints)),paste("TP =",sum(realPpoints))),border="white", col=c(col3,col1),pch=c(1,2),cex=1.5)
legend("topright", legend=c("FP PMV <= 90%","TP PMV <= 90%","FP PMV > 90%","TP PMV > 90%"),border="white", col=c(c(col3,col1),c(col3,col1)),pch=c(c(1,2),c(16,17)),cex=1.2)

resultsub <- sim.sort.result1[1:TOPN,]
sim.logfc <- resultsub[,3]-resultsub[,4]
shapint<-ifelse(resultsub[,1]==1,1,0)
realNpoints <- shapint == 0
realPpoints <- shapint == 1
zeroint <- resultsub[,10]>0.9
plot(1, type="n", xlab="Log Fold Change", ylab="Variance", xlim=c(-4,4), ylim=c(0,4),main="b",cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(sim.logfc[realPpoints], resultsub[,7][realPpoints],col=col1,pch=2,cex =1)
points(sim.logfc[realNpoints], resultsub[,7][realNpoints],col=col3,pch=1)
points(sim.logfc[realPpoints&zeroint], resultsub[,7][realPpoints&zeroint],col=col1,pch=17,cex =1)
points(sim.logfc[realNpoints&zeroint], resultsub[,7][realNpoints&zeroint],col=col3,pch=16)

################################Figure S3#########################################
par(mfrow=c(2,2))
par(oma=c(2,2,0,0))
par(mar=c(3,3,2,1))

pind <- rawpall>0.9
pind2 <- abs(post.var-sim.var)<0.5
pind3 <- abs(t.var-sim.var)<0.6
pind4 <- abs(post.var-sim.var)<0.2
pind5 <- abs(t.var-sim.var)<0.2
cexsize=0.6
col1=rgb(0.1,0.3,0.5,0.8)
col2=rgb(0.8,0.3,0,0.1)
col3=rgb(0.8,0.3,0,0.6)
col4=rgb(0.8,0.3,0,0.3)
col5=rgb(0.5,0.3,0,0.1)
par(mfrow=c(2,2))
par(oma=c(2,2,0,0))
par(mar=c(3,3,2,1))
plot(t.var[pind==FALSE],sim.var[pind==FALSE],xlim= c(0,22),xlab="", ylab="",main = "a",cex.axis=1.4,cex.main=2,pch=16,col=col1,cex =cexsize)
points(t.var[pind==T],sim.var[pind==T],pch=16,col=col2,cex =cexsize)
points(t.var[pind==T&pind3==FALSE],sim.var[pind==T&pind3==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
plot(post.var[pind==FALSE],sim.var[pind==FALSE],xlim= c(0,22),xlab="", ylab="",main = "b",cex.axis=1.4,cex.lab=2,cex.main=2,pch=16,col=col1,cex =cexsize)
points(post.var[pind==T],sim.var[pind==T],pch=16,col=col5,cex =cexsize)
points(post.var[pind==T&pind2==FALSE],sim.var[pind==T&pind2==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
plot(t.var[pind==FALSE],sim.var[pind==FALSE],xlim=c(0,2), ylim=c(0,2),xlab="", ylab="",main = "c",cex.axis=1.4,cex.lab=2,cex.main=2,pch=16,col=col1,cex =cexsize)
points(t.var[pind==T],sim.var[pind==T],pch=16,col=col4,cex =cexsize)
points(t.var[pind==T&pind5==FALSE],sim.var[pind==T&pind5==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
plot(post.var[pind==FALSE],sim.var[pind==FALSE],xlim=c(0,2), ylim=c(0,2),xlab="", ylab="",main = "d",cex.axis=1.4,cex.main=2,pch=16,col=col1,cex =cexsize)
points(post.var[pind==T],sim.var[pind==T],pch=16,col=col4,cex =cexsize)
points(post.var[pind==T&pind4==FALSE],sim.var[pind==T&pind4==FALSE],pch=16,col=col3,cex =cexsize)
abline(0, 1, col="red")
legend("bottomright", legend=c(paste("PMV > 90%"),paste("PMV <= 90%")),border="white", col=c(col3,col1),pch=c(16,16),cex=1.3)

mtext("True Variance", side=2, line=0, cex=2,outer=TRUE)
mtext("Estimated Variance", side=1, line=0, cex=2,outer=TRUE)

################################Figure S4#########################################
par(mfrow=c(2,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))
cexsize=0.6
col1=rgb(0.1,0.3,0.5,0.8)
col2=rgb(0.8,0.3,0,0.1)
col3=rgb(0.8,0.3,0,0.6)
col4=rgb(0.8,0.3,0,0.3)
col5=rgb(0.5,0.3,0,0.1)
par(mfrow=c(2,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))
pind <- rawpall>0.9
pind2<- abs(t.mucontrol-sim.mucontrol)>0.5
pind3<- abs(post.mucontrol.f-sim.mucontrol)>0.5
pind4<- abs(t.pcontrol-sim.pcontrol)>0.1
pind5<- abs(post.pcontrol.f-sim.pcontrol)>0.1
plot(1, type="n", xlab="Estimated Mean", ylab="True Mean",main = "a",xlim=c(0,10), ylim=c(0,10),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(t.mucontrol[pind==FALSE], sim.mucontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(t.mucontrol[pind], sim.mucontrol[pind],col=col2,pch=16,cex=cexsize)
points(t.mucontrol[pind&pind2], sim.mucontrol[pind&pind2],col=col3,pch=16,cex=cexsize)
plot(1, type="n", xlab="Estimated Mean", ylab="True Mean",main = "b",xlim=c(0,10), ylim=c(0,10),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(post.mucontrol.f[pind==FALSE], sim.mucontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(post.mucontrol.f[pind], sim.mucontrol[pind],col=col2,pch=16,cex=cexsize)
points(post.mucontrol.f[pind&pind3], sim.mucontrol[pind&pind3],col=col3,pch=16,cex=cexsize)
#legend("bottomright", legend=c(paste("PMV > 90%"),paste("PMV <= 90%")),border="white", col=c(col3,col1),pch=c(16,16),cex=1.3)

plot(1, type="n", xlab="Estimated BPMV Proprotion", ylab="True BPMV Proportion",main = "c",xlim=c(0,1), ylim=c(0,1),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(t.pcontrol[pind==FALSE], sim.pcontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(t.pcontrol[pind], sim.pcontrol[pind],col=col2,pch=16,cex=cexsize)
points(t.pcontrol[pind&pind4], sim.pcontrol[pind&pind4],col=col3,pch=16,cex=cexsize)
plot(1, type="n", xlab="Estimated BPMV Proprotion", ylab="True BPMV Proportion",main = "d",xlim=c(0,1), ylim=c(0,1),cex.axis=1.4,cex.lab=1.8,cex.main=2)
points(post.pcontrol.f[pind==FALSE], sim.pcontrol[pind==FALSE],col=col1,pch=16,cex=cexsize)
points(post.pcontrol.f[pind], sim.pcontrol[pind],col=col2,pch=16,cex=cexsize)
points(post.pcontrol.f[pind&pind5], sim.pcontrol[pind&pind5],col=col3,pch=16,cex=cexsize)
legend("bottomright", legend=c(paste("PMV > 90%"),paste("PMV <= 90%")),border="white", col=c(col3,col1),pch=c(16,16),cex=1.3)

################################Figure S5#########################################
par(mfrow=c(2,2))
par(oma=c(0,0,0,0))
par(mar=c(5,5,2,1))

numobs=100
DASEV.TRP <-c()
DASEV.TRP.M <-c()
LIM.TRP <-c()
LIM.TRP.M <-c()
DASEV.FDR <-c()
DASEV.FDR.M <-c()
LIM.FDR <-c()
LIM.FDR.M <-c()
a1_all <- c()
b1_all <- c()
c1_all <- c()
d1_all <- c()
a2_all <- c()
b2_all <- c()
c2_all <- c()
d2_all <- c()
a3_all <- c()
b3_all <- c()
c3_all <- c()
d3_all <- c()
for (k in 1:100){
  
  fname <- paste("C:\\DASEV\\S2\\simresult_N_",numobs,"_",k,".Rdata",sep="")
  fname1 <- paste("C:\\DASEV\\S2\\otherresult_N_",numobs,"_",k,".Rdata",sep="")
  fnamep <- paste("C:\\DASEV\\S2\\paradata_N_",numobs,"_",k,".Rdata",sep="")
  fname.data <- paste("C:\\DASEV\\S2\\simdata_N_",numobs,"_",k,".Rdata",sep="")
  
  load(file=fname.data) # data2
  datatmp <- data2
  count <- rowSums(datatmp==0)
  rawpall  <- count/(2*numobs)
  
  load(file=fname) #sim.result
  est <- sim.result$estimates
  p <-  sim.result$pvalue_both
  p.m <-  sim.result$pvalue_mean
  p.p <-  sim.result$pvalue_zero
  post.mucontrol.f <- est[,1]
  post.mucase.f <- est[,1]+est[,2]
  post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
  post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
  post.var <- est[,5]^2
  
  load(file=fname1) #other.result
  t.p <- other.result[,1]
  t.var <- other.result[,2]
  t.mucontrol <- other.result[,3]
  t.mucase <- other.result[,4]
  t.pcontrol <- other.result[,5]
  t.pcase <- other.result[,6]
  t.p.m <- other.result[,7]
  t.p.p <- other.result[,13]
  
  load(file=fnamep) #parameters
  sim.lod <- paradata[,1]
  sim.mucontrol <- paradata[,2]
  sim.pcontrol <- paradata[,3]
  sim.sd <- paradata[,4]
  sim.mucase <- paradata[,5]
  sim.pcase <- paradata[,6]
  sim.geneid <- paradata[,7]
  sim.geneid1 <- paradata[,8]
  numint <- paradata[,9]
  
  sim.var <- sim.sd^2
  
  q<-p.adjust(p,method='BH')
  qmean<-p.adjust(p.m,method='BH')
  qzero<-p.adjust(p.p,method='BH')
  t.q<-p.adjust(t.p,method='BH')
  t.qmean<-p.adjust(t.p.m,method="BH")
  t.qzero<-p.adjust(t.p.p,method='BH')
  
  rank <- c(1:length(sim.geneid))
  
  a1_ <- sum(qzero<=0.01 & sim.geneid==1)
  b1_ <- sum(qzero<=0.01 & sim.geneid==0)
  c1_ <- sum(t.qzero<=0.01 & sim.geneid==1)
  d1_ <- sum(t.qzero<=0.01 & sim.geneid==0)
  
  a2_ <- sum(qzero<=0.05 & sim.geneid==1)
  b2_ <- sum(qzero<=0.05 & sim.geneid==0)
  c2_ <- sum(t.qzero<=0.05 & sim.geneid==1)
  d2_ <- sum(t.qzero<=0.05 & sim.geneid==0)
  
  a3_ <- sum(qzero<=0.1 & sim.geneid==1)
  b3_ <- sum(qzero<=0.1 & sim.geneid==0)
  c3_ <- sum(t.qzero<=0.1 & sim.geneid==1)
  d3_ <- sum(t.qzero<=0.1 & sim.geneid==0)
  
  a1_all <- c(a1_all,a1_)
  b1_all <- c(b1_all,b1_)
  c1_all <- c(c1_all,c1_)
  d1_all <- c(d1_all,d1_)
  a2_all <- c(a2_all,a2_)
  b2_all <- c(b2_all,b2_)
  c2_all <- c(c2_all,c2_)
  d2_all <- c(d2_all,d2_)
  a3_all <- c(a3_all,a3_)
  b3_all <- c(b3_all,b3_)
  c3_all <- c(c3_all,c3_)
  d3_all <- c(d3_all,d3_)
  
  sim.sort.result <- cbind(sim.geneid,q,post.mucase.f,numint,p)
  sim.sort.result <-sim.sort.result[order(p),]
  sim.sort.result <- cbind(sim.sort.result,rank)
  sim.sort.result1 <- cbind(sim.geneid,qzero,post.mucase.f,numint,p.p)
  sim.sort.result1 <-sim.sort.result1[order(p.p),]
  sim.sort.result1 <- cbind(sim.sort.result1,rank)
  t.sort.result <- cbind(sim.geneid,t.q,t.mucase,numint,t.p)
  t.sort.result <-t.sort.result[order(t.p),]
  t.sort.result <- cbind(t.sort.result,rank)
  t.sort.result1 <- cbind(sim.geneid,t.qzero,t.mucase,numint,t.p.p)
  t.sort.result1 <-t.sort.result1[order(t.p.p),]
  t.sort.result1 <- cbind(t.sort.result1,rank)
  
  #tOP RANKED FEATURES AND TRUE POSITIVE RATE
  TOPn <- seq(5,2000,5)
  TPR <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR <- c(TPR,tpr)
    TPR[is.na(TPR)] <- 0
  }
  DASEV.TRP <-cbind(DASEV.TRP,TPR)
  
  TPR.M <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR.M <- c(TPR.M,tpr)
    TPR.M[is.na(TPR.M)] <- 0
  }
  DASEV.TRP.M <-cbind(DASEV.TRP.M,TPR.M)
  
  T.TPR <- c()
  for (i in TOPn){
    resultsub <- t.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR <- c(T.TPR,tpr)
    T.TPR[is.na(T.TPR)] <- 0
  }
  LIM.TRP <- cbind(LIM.TRP,T.TPR)
  
  T.TPR.M <- c()
  for (i in TOPn){
    resultsub <- t.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR.M <- c(T.TPR.M,tpr)
    T.TPR.M[is.na(T.TPR.M)] <- 0
  }
  LIM.TRP.M <- cbind(LIM.TRP.M,T.TPR.M)
  
  
  cutoffs <- seq(0,1,0.005)
  
  FDR <- c()
  FDR.M <- c()
  T.FDR <- c()
  T.FDR.M <- c()
  for (j in cutoffs){
    DE=sum(q<j & sim.geneid==1)
    TOT=sum(q<j)
    if (TOT==0) fdr=0 else fdr = 1- DE/TOT
    FDR <- c(FDR,fdr)
    
    DE.M=sum(qmean<j & sim.geneid==1)
    TOT.M=sum(qmean<j)
    if (TOT.M==0) fdr=0 else fdr = 1- DE.M/TOT.M
    FDR.M <- c(FDR.M,fdr)
    
    T.DE=sum(t.q<j & sim.geneid==1)
    T.TOT=sum(t.q<j)
    if (T.TOT==0) fdr=0 else fdr = 1- T.DE/T.TOT
    T.FDR <- c(T.FDR,fdr)
    
    T.DE.M=sum(t.qmean<j & sim.geneid==1)
    T.TOT.M=sum(t.qmean<j)
    if (T.TOT.M==0) fdr=0 else fdr = 1- T.DE.M/T.TOT.M
    T.FDR.M <- c(T.FDR.M,fdr)
  }
  DASEV.FDR <-cbind(DASEV.FDR,FDR)
  DASEV.FDR.M <-cbind(DASEV.FDR.M,FDR.M)
  LIM.FDR <-cbind(LIM.FDR,T.FDR)
  LIM.FDR.M <-cbind(LIM.FDR.M,T.FDR.M)
}

mean.DASEV.TRP_100 <-rowMeans(DASEV.TRP)
mean.DASEV.TRP.M_100  <-rowMeans(DASEV.TRP.M)
mean.LIM.TRP_100  <-rowMeans(LIM.TRP)
mean.LIM.TRP.M_100  <-rowMeans(LIM.TRP.M)
mean.DASEV.FDR_100  <-rowMeans(DASEV.FDR)
mean.DASEV.FDR.M_100  <-rowMeans(DASEV.FDR.M)
mean.LIM.FDR_100  <-rowMeans(LIM.FDR)
mean.LIM.FDR.M_100  <-rowMeans(LIM.FDR.M)
a1<- round(mean(a1_all))
b1<- round(mean(b1_all))
c1<- round(mean(c1_all))
d1<- round(mean(d1_all))
a2<- round(mean(a2_all))
b2<- round(mean(b2_all))
c2<- round(mean(c2_all))
d2<- round(mean(d2_all))
a3<- round(mean(a3_all))
b3<- round(mean(b3_all))
c3<- round(mean(c3_all))
d3<- round(mean(d3_all))

v1_100  <- c(a1,a2,a3)
v2_100  <- c(-b1,-b2,-b3)
v3_100  <- c(c1,c2,c3)
v4_100  <- c(-d1,-d2,-d3)

tfdrd1 <- b1_all/(b1_all+a1_all)
tfdrd2 <- b2_all/(b2_all+a2_all)
tfdrd3 <- b3_all/(b3_all+a3_all)
tfdrt1 <- d1_all/(d1_all+c1_all)
tfdrt2 <- d2_all/(d2_all+c2_all)
tfdrt3 <- d3_all/(d3_all+c3_all)

tprrate_100_2  <- c(mean(tfdrd1[!is.na(tfdrd1)]),mean(tfdrt1[!is.na(tfdrt1)]),
                    mean(tfdrd2[!is.na(tfdrd2)]),mean(tfdrt2[!is.na(tfdrt2)]),
                    mean(tfdrd3[!is.na(tfdrd3)]),mean(tfdrt3[!is.na(tfdrd3)]))
rate_100 <- paste(round(100*tprrate_100_2, 2), "%", sep="")


numobs=200
DASEV.TRP <-c()
DASEV.TRP.M <-c()
LIM.TRP <-c()
LIM.TRP.M <-c()
DASEV.FDR <-c()
DASEV.FDR.M <-c()
LIM.FDR <-c()
LIM.FDR.M <-c()
a1_all <- c()
b1_all <- c()
c1_all <- c()
d1_all <- c()
a2_all <- c()
b2_all <- c()
c2_all <- c()
d2_all <- c()
a3_all <- c()
b3_all <- c()
c3_all <- c()
d3_all <- c()
for (k in 1:100){
  fname <- paste("C:\\DASEV\\S2\\simresult_N_",numobs,"_",k,".Rdata",sep="")
  fname1 <- paste("C:\\DASEV\\S2\\otherresult_N_",numobs,"_",k,".Rdata",sep="")
  fnamep <- paste("C:\\DASEV\\S2\\paradata_N_",numobs,"_",k,".Rdata",sep="")
  fname.data <- paste("C:\\DASEV\\S2\\simdata_N_",numobs,"_",k,".Rdata",sep="")
  
  load(file=fname.data) # data2
  datatmp <- data2
  count <- rowSums(datatmp==0)
  rawpall  <- count/(2*numobs)
  
  load(file=fname) #sim.result
  est <- sim.result$estimates
  p <-  sim.result$pvalue_both
  p.m <-  sim.result$pvalue_mean
  p.p <-  sim.result$pvalue_zero
  post.mucontrol.f <- est[,1]
  post.mucase.f <- est[,1]+est[,2]
  post.pcontrol.f <- 1 - exp(est[,3])/(1+exp(est[,3]))
  post.pcase.f <- 1 - exp(est[,3]+ est[,4])/(1+exp(est[,3])+est[,4])
  post.var <- est[,5]^2
  
  load(file=fname1) #other.result
  t.p <- other.result[,1]
  t.var <- other.result[,2]
  t.mucontrol <- other.result[,3]
  t.mucase <- other.result[,4]
  t.pcontrol <- other.result[,5]
  t.pcase <- other.result[,6]
  t.p.m <- other.result[,7]
  t.p.p <- other.result[,13]
  
  load(file=fnamep) #parameters
  sim.lod <- paradata[,1]
  sim.mucontrol <- paradata[,2]
  sim.pcontrol <- paradata[,3]
  sim.sd <- paradata[,4]
  sim.mucase <- paradata[,5]
  sim.pcase <- paradata[,6]
  sim.geneid <- paradata[,7]
  sim.geneid1 <- paradata[,8]
  numint <- paradata[,9]
  
  sim.var <- sim.sd^2
  
  q<-p.adjust(p,method='BH')
  qmean<-p.adjust(p.m,method='BH')
  qzero<-p.adjust(p.p,method='BH')
  t.q<-p.adjust(t.p,method='BH')
  t.qmean<-p.adjust(t.p.m,method="BH")
  t.qzero<-p.adjust(t.p.p,method='BH')
  
  rank <- c(1:length(sim.geneid))
  
  
  a1_ <- sum(qzero<=0.01 & sim.geneid==1)
  b1_ <- sum(qzero<=0.01 & sim.geneid==0)
  c1_ <- sum(t.qzero<=0.01 & sim.geneid==1)
  d1_ <- sum(t.qzero<=0.01 & sim.geneid==0)
  
  a2_ <- sum(qzero<=0.05 & sim.geneid==1)
  b2_ <- sum(qzero<=0.05 & sim.geneid==0)
  c2_ <- sum(t.qzero<=0.05 & sim.geneid==1)
  d2_ <- sum(t.qzero<=0.05 & sim.geneid==0)
  
  a3_ <- sum(qzero<=0.1 & sim.geneid==1)
  b3_ <- sum(qzero<=0.1 & sim.geneid==0)
  c3_ <- sum(t.qzero<=0.1 & sim.geneid==1)
  d3_ <- sum(t.qzero<=0.1 & sim.geneid==0)
  
  a1_all <- c(a1_all,a1_)
  b1_all <- c(b1_all,b1_)
  c1_all <- c(c1_all,c1_)
  d1_all <- c(d1_all,d1_)
  a2_all <- c(a2_all,a2_)
  b2_all <- c(b2_all,b2_)
  c2_all <- c(c2_all,c2_)
  d2_all <- c(d2_all,d2_)
  a3_all <- c(a3_all,a3_)
  b3_all <- c(b3_all,b3_)
  c3_all <- c(c3_all,c3_)
  d3_all <- c(d3_all,d3_)
  
  sim.sort.result <- cbind(sim.geneid,q,post.mucase.f,numint,p)
  sim.sort.result <-sim.sort.result[order(p),]
  sim.sort.result <- cbind(sim.sort.result,rank)
  sim.sort.result1 <- cbind(sim.geneid,qzero,post.mucase.f,numint,p.p)
  sim.sort.result1 <-sim.sort.result1[order(p.p),]
  sim.sort.result1 <- cbind(sim.sort.result1,rank)
  t.sort.result <- cbind(sim.geneid,t.q,t.mucase,numint,t.p)
  t.sort.result <-t.sort.result[order(t.p),]
  t.sort.result <- cbind(t.sort.result,rank)
  t.sort.result1 <- cbind(sim.geneid,t.qzero,t.mucase,numint,t.p.p)
  t.sort.result1 <-t.sort.result1[order(t.p.p),]
  t.sort.result1 <- cbind(t.sort.result1,rank)
  
  #tOP RANKED FEATURES AND TRUE POSITIVE RATE
  TOPn <- seq(5,2000,5)
  TPR <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR <- c(TPR,tpr)
    TPR[is.na(TPR)] <- 0
  }
  DASEV.TRP <-cbind(DASEV.TRP,TPR)
  
  TPR.M <- c()
  for (i in TOPn){
    resultsub <- sim.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    TPR.M <- c(TPR.M,tpr)
    TPR.M[is.na(TPR.M)] <- 0
  }
  DASEV.TRP.M <-cbind(DASEV.TRP.M,TPR.M)
  
  T.TPR <- c()
  for (i in TOPn){
    resultsub <- t.sort.result[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR <- c(T.TPR,tpr)
    T.TPR[is.na(T.TPR)] <- 0
  }
  LIM.TRP <- cbind(LIM.TRP,T.TPR)
  
  T.TPR.M <- c()
  for (i in TOPn){
    resultsub <- t.sort.result1[1:i,]
    DE <- resultsub[,1][resultsub[,1]==1]
    DE[is.na(DE)] <- 0
    tpr <- length(DE)/nrow(resultsub)
    T.TPR.M <- c(T.TPR.M,tpr)
    T.TPR.M[is.na(T.TPR.M)] <- 0
  }
  LIM.TRP.M <- cbind(LIM.TRP.M,T.TPR.M)
  
  
  cutoffs <- seq(0,1,0.005)
  
  FDR <- c()
  FDR.M <- c()
  T.FDR <- c()
  T.FDR.M <- c()
  for (j in cutoffs){
    DE=sum(q<j & sim.geneid==1)
    TOT=sum(q<j)
    if (TOT==0) fdr=0 else fdr = 1- DE/TOT
    FDR <- c(FDR,fdr)
    
    DE.M=sum(qmean<j & sim.geneid==1)
    TOT.M=sum(qmean<j)
    if (TOT.M==0) fdr=0 else fdr = 1- DE.M/TOT.M
    FDR.M <- c(FDR.M,fdr)
    
    T.DE=sum(t.q<j & sim.geneid==1)
    T.TOT=sum(t.q<j)
    if (T.TOT==0) fdr=0 else fdr = 1- T.DE/T.TOT
    T.FDR <- c(T.FDR,fdr)
    
    T.DE.M=sum(t.qmean<j & sim.geneid==1)
    T.TOT.M=sum(t.qmean<j)
    if (T.TOT.M==0) fdr=0 else fdr = 1- T.DE.M/T.TOT.M
    T.FDR.M <- c(T.FDR.M,fdr)
  }
  DASEV.FDR <-cbind(DASEV.FDR,FDR)
  DASEV.FDR.M <-cbind(DASEV.FDR.M,FDR.M)
  LIM.FDR <-cbind(LIM.FDR,T.FDR)
  LIM.FDR.M <-cbind(LIM.FDR.M,T.FDR.M)
}
mean.DASEV.TRP <-rowMeans(DASEV.TRP)
mean.DASEV.TRP.M <-rowMeans(DASEV.TRP.M)
mean.LIM.TRP <-rowMeans(LIM.TRP)
mean.LIM.TRP.M <-rowMeans(LIM.TRP.M)
mean.DASEV.FDR <-rowMeans(DASEV.FDR)
mean.DASEV.FDR.M <-rowMeans(DASEV.FDR.M)
mean.LIM.FDR <-rowMeans(LIM.FDR)
mean.LIM.FDR.M <-rowMeans(LIM.FDR.M)
a1<- round(mean(a1_all))
b1<- round(mean(b1_all))
c1<- round(mean(c1_all))
d1<- round(mean(d1_all))
a2<- round(mean(a2_all))
b2<- round(mean(b2_all))
c2<- round(mean(c2_all))
d2<- round(mean(d2_all))
a3<- round(mean(a3_all))
b3<- round(mean(b3_all))
c3<- round(mean(c3_all))
d3<- round(mean(d3_all))

plot(TOPn,mean.DASEV.TRP.M_100 ,type="l",col="darkorange2", xlim=c(0,300),ylim=c(0.4,1), xlab="Top Ranked Features", ylab="True Positive Rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="a",lwd=2)
lines(TOPn,mean.LIM.TRP.M_100 ,type="l",col="royalblue4", lwd=2)
legend("bottomleft", legend=c("DASEV","TLK"),border="white", col=c('darkorange2', 'royalblue4'),lty=c(1,1),cex=1.2)

plot(TOPn,mean.DASEV.TRP.M,type="l",col="darkorange2", xlim=c(0,300),ylim=c(0.4,1), xlab="Top Ranked Features", ylab="True Positive Rate",cex.axis=1.4,cex.lab=1.8,cex.main=2,main ="b",lwd=2)
lines(TOPn,mean.LIM.TRP.M,type="l",col="royalblue4", lwd=2)
legend("bottomleft", legend=c("DASEV","TLK"),border="white", col=c('darkorange2', 'royalblue4'),lty=c(1,1),cex=1.2)

barplot(matrix(c(v1_100 ,v3_100 ),nrow=2, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4'),ylim=c(-50,50),xlab="Reported FDR", ylab = "Counts",main = "c",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(v2_100 ,v4_100 ),nrow=2, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2'), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,4.5,5.5,7.5,8.5), y=(c(20,25,30,35,40,45)-15), labels=rate_100 ,cex=1.2)
legend("bottomleft", legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2"),cex=1.2,pch=c(15,15,15,15))

v1 <- c(a1,a2,a3)
v2 <- c(-b1,-b2,-b3)
v3 <- c(c1,c2,c3)
v4 <- c(-d1,-d2,-d3)
cutoffs2<- c(0.01,0.05,0.1)
tfdrd1 <- b1_all/(b1_all+a1_all)
tfdrd2 <- b2_all/(b2_all+a2_all)
tfdrd3 <- b3_all/(b3_all+a3_all)
tfdrt1 <- d1_all/(d1_all+c1_all)
tfdrt2 <- d2_all/(d2_all+c2_all)
tfdrt3 <- d3_all/(d3_all+c3_all)

tprrate_2  <- c(mean(tfdrd1[!is.na(tfdrd1)]),mean(tfdrt1[!is.na(tfdrt1)]),
                mean(tfdrd2[!is.na(tfdrd2)]),mean(tfdrt2[!is.na(tfdrt2)]),
                mean(tfdrd3[!is.na(tfdrd3)]),mean(tfdrt3[!is.na(tfdrd3)]))
rate <- paste(round(100*tprrate_2, 2), "%", sep="")
barplot(matrix(c(v1,v3),nrow=2, byrow=T),beside=T, names.arg =c("1%","5%","10%"),col=c('darkorange2', 'royalblue4'),ylim=c(-100,150),xlab="Reported FDR", ylab = "Counts",main = "d",cex.axis=1.4,cex.names=1.4,cex.lab=1.8,cex.main=2)
barplot(matrix(c(v2,v4),nrow=2, byrow=T),beside=T, col=c('lightgoldenrod2','dodgerblue2'), add=T,cex.axis=1.4)
text(x=c(1.5,2.5,4.5,5.5,7.5,8.5), y=(c(150,165,200,215,230,245)-100), labels=rate,cex=1.2)
legend("bottomleft", legend=c("DASEV TP","DASEV FP","TLK TP","TLK FP"),border="white", col=c("darkorange2","lightgoldenrod2","royalblue4","dodgerblue2"),cex=1.2,pch=c(15,15,15,15))